Load packages

library(ggplot2)
library(cowplot)
library(lme4)
library(pbkrtest)
library(tidyr)
library(dplyr)
library(optimx)
library(broom)

Data description

Seven censuses were conducted in this experiment (from 14fa to 17fa, twice a year). Only those seedlings inside exclosures were used in the analyses. Treatment: C (no manipulation), S (snow removal treatment).

Initial analysis

getwd()
## [1] "E:/R/My code/Git/Changbaishan/Changbaishan1"
## Read in data
snowdat.ag_adult <- read.csv("data/snowdat.ag_adult.csv")

snowdat.ag_adult$census <- factor(snowdat.ag_adult$census)
snowdat.ag_adult$site <- factor(snowdat.ag_adult$site)
snowdat.ag_adult$census <- factor(snowdat.ag_adult$census, levels=c('15sp', '15fa', '16sp', '16fa','17sp','17fa'))
summary(snowdat.ag_adult)
##        X          site       quadrat         sp.      treatment
##  Min.   :   1.0   1:445   Min.   :102   FRMA   :217   C:657    
##  1st Qu.: 291.2   2:455   1st Qu.:203   TIAM   :189   S:505    
##  Median : 581.5   3:262   Median :306   PHSC   :115            
##  Mean   : 581.5           Mean   :299   EUMA   : 82            
##  3rd Qu.: 871.8           3rd Qu.:404   DEGL   : 60            
##  Max.   :1162.0           Max.   :508   ACPS   : 52            
##                                         (Other):447            
##  growth.form  census        deaths            total       
##  shrub:467   15sp:208   Min.   : 0.0000   Min.   : 1.000  
##  tree :695   15fa:117   1st Qu.: 0.0000   1st Qu.: 1.000  
##              16sp: 94   Median : 0.0000   Median : 2.000  
##              16fa:225   Mean   : 0.6945   Mean   : 2.543  
##              17sp:219   3rd Qu.: 1.0000   3rd Qu.: 3.000  
##              17fa:299   Max.   :17.0000   Max.   :21.000  
##                                                           
##      survs         quad.unique          quad.sp         A.sum     
##  Min.   : 0.000   1_1_308: 43   1_1_102_EUAL:   6   Min.   :4487  
##  1st Qu.: 1.000   2_1_403: 43   1_1_102_PHSC:   6   1st Qu.:5335  
##  Median : 1.000   2_1_501: 43   1_1_104_DEGL:   6   Median :5708  
##  Mean   : 1.849   2_1_104: 37   1_1_104_EUMA:   6   Mean   :5902  
##  3rd Qu.: 2.000   1_1_507: 35   1_1_203_DEGL:   6   3rd Qu.:6533  
##  Max.   :14.000   1_1_409: 32   1_1_205_PHSC:   6   Max.   :9223  
##                   (Other):929   (Other)     :1126                 
##      A.con              A.het         con.dens        A.con_curt    
##  Min.   :   0.000   Min.   :2827   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:   0.000   1st Qu.:4766   1st Qu.: 0.000   1st Qu.: 0.000  
##  Median :   3.047   Median :5455   Median : 1.000   Median : 1.450  
##  Mean   : 446.019   Mean   :5456   Mean   : 4.998   Mean   : 4.357  
##  3rd Qu.: 648.325   3rd Qu.:6120   3rd Qu.: 8.000   3rd Qu.: 8.655  
##  Max.   :4003.948   Max.   :9223   Max.   :41.000   Max.   :15.879  
##                                                                     
##    A.het_curt      A.sum_curt    con.dens_curt  
##  Min.   :14.14   Min.   :16.49   Min.   :0.000  
##  1st Qu.:16.83   1st Qu.:17.47   1st Qu.:0.000  
##  Median :17.60   Median :17.87   Median :1.000  
##  Mean   :17.53   Mean   :18.03   Mean   :1.010  
##  3rd Qu.:18.29   3rd Qu.:18.69   3rd Qu.:2.000  
##  Max.   :20.97   Max.   :20.97   Max.   :3.448  
## 
###------------------------ Initial plots -------------------------###
ggplot(data=snowdat.ag_adult, 
       aes(x=treatment, y=survs/total, colour=as.factor(census))) + 
  geom_boxplot() +
  facet_wrap(~growth.form)

group_by(snowdat.ag_adult, treatment, growth.form, census) %>% 
  summarise(surv = mean(survs/total), 
            se.survs=sd(survs/total)/sqrt(length(survs/total))) %>% 
  ggplot(aes(x=treatment, y=surv, ymin=surv-se.survs, ymax=surv+se.survs, 
             colour=as.factor(census))) + geom_errorbar() + geom_point() +
  facet_wrap(~growth.form) 

group_by(snowdat.ag_adult, treatment, growth.form, census) %>% 
  summarise(surv = mean(survs/total), 
            se.survs=sd(survs/total)/sqrt(length(survs/total))) %>% 
  ggplot(aes(x=census, y=surv, ymin=surv-se.survs, ymax=surv+se.survs, 
             colour=as.factor(treatment))) + geom_errorbar() + geom_point() +
  facet_wrap(~growth.form) 

###----------------------- Survival models ------------------------###
snowdat.ag_adult$site <- factor(snowdat.ag_adult$site)
contrasts(snowdat.ag_adult$site) <- "contr.sum"

#### tree
## with or without neighboring trees
mod.snow.tree.surv1 <- glmer(cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
                               (1|quad.unique) + (1|sp.) + (1|census),
                             data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
# without neighboring trees
mod.snow.tree.surv2 <- glmer(cbind(survs, deaths) ~ site + treatment +
                               (1|quad.unique) + (1|sp.) + (1|census),
                             data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.tree.surv1, mod.snow.tree.surv2)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv2: cbind(survs, deaths) ~ site + treatment + (1 | quad.unique) + 
## mod.snow.tree.surv2:     (1 | sp.) + (1 | census)
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt + 
## mod.snow.tree.surv1:     (1 | quad.unique) + (1 | sp.) + (1 | census)
##                     Df    AIC  BIC  logLik deviance  Chisq Chi Df
## mod.snow.tree.surv2  7 1343.2 1375 -664.62   1329.2              
## mod.snow.tree.surv1  9 1344.1 1385 -663.07   1326.1 3.0972      2
##                     Pr(>Chisq)
## mod.snow.tree.surv2           
## mod.snow.tree.surv1     0.2125
# adding area of neighboring trees make no sense to the model, but I would like to keep it as the same as the survival analysis of the pesticide experiment.

# model sensitivity tests
mod.snow.tree.surv.diags <- augment(mod.snow.tree.surv2)
qplot(data=mod.snow.tree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data=mod.snow.tree.surv.diags, x=treatment, y=.wtres, geom="boxplot") 

library(car)
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
qqPlot(resid(mod.snow.tree.surv1))

library(sjPlot)
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
sjp.lmer(mod.snow.tree.surv2, type='fe')

sjp.lmer(mod.snow.tree.surv2, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.snow.tree.surv2, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment + (1 | quad.unique) +  
##     (1 | sp.) + (1 | census)
##    Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1343.2   1375.0   -664.6   1329.2      688 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2984 -0.6223  0.2682  0.6565  3.0701 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1182   0.3437  
##  sp.         (Intercept) 0.7354   0.8576  
##  census      (Intercept) 1.5785   1.2564  
## Number of obs: 695, groups:  quad.unique, 60; sp., 15; census, 6
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.38884    0.58816   0.661    0.509
## site1        0.14549    0.10515   1.384    0.166
## site2       -0.04758    0.10295  -0.462    0.644
## treatmentS  -0.07413    0.14617  -0.507    0.612
##-----------------------------##
# Snow removal treatment didn't change the survival of seedlings, 
# but the BLUP values of census show that snow removal duration 
# may has a potential relationship with survival.

#### duration
## duration is presented as months from the beginning of the experiment
snowdat.ag_adult$duration <- ifelse(snowdat.ag_adult$census=="15sp", 9,
                                    ifelse(snowdat.ag_adult$census=="15fa", 12,
                                           ifelse(snowdat.ag_adult$census=="16sp", 21,
                                                  ifelse(snowdat.ag_adult$census=="16fa", 24,
                                                         ifelse(snowdat.ag_adult$census=="17sp", 32,
                                                                ifelse(snowdat.ag_adult$census=="17fa", 35, NA))))))
summary(snowdat.ag_adult)
##        X          site       quadrat         sp.      treatment
##  Min.   :   1.0   1:445   Min.   :102   FRMA   :217   C:657    
##  1st Qu.: 291.2   2:455   1st Qu.:203   TIAM   :189   S:505    
##  Median : 581.5   3:262   Median :306   PHSC   :115            
##  Mean   : 581.5           Mean   :299   EUMA   : 82            
##  3rd Qu.: 871.8           3rd Qu.:404   DEGL   : 60            
##  Max.   :1162.0           Max.   :508   ACPS   : 52            
##                                         (Other):447            
##  growth.form  census        deaths            total       
##  shrub:467   15sp:208   Min.   : 0.0000   Min.   : 1.000  
##  tree :695   15fa:117   1st Qu.: 0.0000   1st Qu.: 1.000  
##              16sp: 94   Median : 0.0000   Median : 2.000  
##              16fa:225   Mean   : 0.6945   Mean   : 2.543  
##              17sp:219   3rd Qu.: 1.0000   3rd Qu.: 3.000  
##              17fa:299   Max.   :17.0000   Max.   :21.000  
##                                                           
##      survs         quad.unique          quad.sp         A.sum     
##  Min.   : 0.000   1_1_308: 43   1_1_102_EUAL:   6   Min.   :4487  
##  1st Qu.: 1.000   2_1_403: 43   1_1_102_PHSC:   6   1st Qu.:5335  
##  Median : 1.000   2_1_501: 43   1_1_104_DEGL:   6   Median :5708  
##  Mean   : 1.849   2_1_104: 37   1_1_104_EUMA:   6   Mean   :5902  
##  3rd Qu.: 2.000   1_1_507: 35   1_1_203_DEGL:   6   3rd Qu.:6533  
##  Max.   :14.000   1_1_409: 32   1_1_205_PHSC:   6   Max.   :9223  
##                   (Other):929   (Other)     :1126                 
##      A.con              A.het         con.dens        A.con_curt    
##  Min.   :   0.000   Min.   :2827   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:   0.000   1st Qu.:4766   1st Qu.: 0.000   1st Qu.: 0.000  
##  Median :   3.047   Median :5455   Median : 1.000   Median : 1.450  
##  Mean   : 446.019   Mean   :5456   Mean   : 4.998   Mean   : 4.357  
##  3rd Qu.: 648.325   3rd Qu.:6120   3rd Qu.: 8.000   3rd Qu.: 8.655  
##  Max.   :4003.948   Max.   :9223   Max.   :41.000   Max.   :15.879  
##                                                                     
##    A.het_curt      A.sum_curt    con.dens_curt      duration   
##  Min.   :14.14   Min.   :16.49   Min.   :0.000   Min.   : 9.0  
##  1st Qu.:16.83   1st Qu.:17.47   1st Qu.:0.000   1st Qu.:12.0  
##  Median :17.60   Median :17.87   Median :1.000   Median :24.0  
##  Mean   :17.53   Mean   :18.03   Mean   :1.010   Mean   :24.2  
##  3rd Qu.:18.29   3rd Qu.:18.69   3rd Qu.:2.000   3rd Qu.:35.0  
##  Max.   :20.97   Max.   :20.97   Max.   :3.448   Max.   :35.0  
## 
#### tree
mod.snow.tree.surv3 <- glmer(cbind(survs, deaths) ~ site + treatment*duration +
                               A.con_curt + A.het_curt +
                               (1|quad.unique) + (1|sp.),
                             data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.snow.tree.surv4 <- glmer(cbind(survs, deaths) ~ treatment*duration +
                               A.con_curt + A.het_curt +
                               (1|quad.unique) + (1|census) + (1|sp.),
                             data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.snow.tree.surv1, mod.snow.tree.surv3, mod.snow.tree.surv4)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt + 
## mod.snow.tree.surv1:     (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.tree.surv4: cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt + 
## mod.snow.tree.surv4:     (1 | quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.surv3: cbind(survs, deaths) ~ site + treatment * duration + A.con_curt + 
## mod.snow.tree.surv3:     A.het_curt + (1 | quad.unique) + (1 | sp.)
##                     Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.snow.tree.surv1  9 1344.1 1385.0 -663.07   1326.1              
## mod.snow.tree.surv4  9 1285.9 1326.8 -633.93   1267.9 58.273      0
## mod.snow.tree.surv3 10 1428.2 1473.6 -704.09   1408.2  0.000      1
##                     Pr(>Chisq)    
## mod.snow.tree.surv1               
## mod.snow.tree.surv4     <2e-16 ***
## mod.snow.tree.surv3          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.tree.surv.diags <- augment(mod.snow.tree.surv4)
qplot(data=mod.snow.tree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

sjp.lmer(mod.snow.tree.surv4, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.snow.tree.surv4))

summary(mod.snow.tree.surv4, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt +  
##     (1 | quad.unique) + (1 | census) + (1 | sp.)
##    Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1285.9   1326.8   -633.9   1267.9      686 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1175 -0.5715  0.2764  0.6530  2.6078 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.09238  0.3039  
##  sp.         (Intercept) 0.49931  0.7066  
##  census      (Intercept) 0.36476  0.6040  
## Number of obs: 695, groups:  quad.unique, 60; sp., 15; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -2.10421    1.40830  -1.494  0.13514    
## treatmentS          -2.86182    0.43677  -6.552 5.67e-11 ***
## duration             0.08102    0.02725   2.973  0.00295 ** 
## A.con_curt          -0.03530    0.02358  -1.497  0.13439    
## A.het_curt           0.06162    0.06669   0.924  0.35552    
## treatmentS:duration  0.10083    0.01482   6.804 1.02e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## only according to AIC values, model 4 is better. But the model looks problemic.Other approaches?


#### season
snowdat.ag_adult$season <- ifelse(snowdat.ag_adult$census=='15fa'| 
                                    snowdat.ag_adult$census=='16fa'|
                                    snowdat.ag_adult$census=='17fa', 'fa', 'sp')
mod.snow.tree.surv5 <- glmer(cbind(survs, deaths) ~ site + treatment*season +
                               A.con_curt + A.het_curt +
                               (1|quad.unique) + (1|census) + (1|sp.),
                             data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                             control=glmerControl(optimizer="optimx", 
                                                  optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.snow.tree.surv1, mod.snow.tree.surv3, mod.snow.tree.surv5)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt + 
## mod.snow.tree.surv1:     (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.tree.surv3: cbind(survs, deaths) ~ site + treatment * duration + A.con_curt + 
## mod.snow.tree.surv3:     A.het_curt + (1 | quad.unique) + (1 | sp.)
## mod.snow.tree.surv5: cbind(survs, deaths) ~ site + treatment * season + A.con_curt + 
## mod.snow.tree.surv5:     A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
##                     Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.snow.tree.surv1  9 1344.1 1385.0 -663.07   1326.1              
## mod.snow.tree.surv3 10 1428.2 1473.6 -704.09   1408.2   0.00      1
## mod.snow.tree.surv5 11 1308.0 1358.0 -643.01   1286.0 122.17      1
##                     Pr(>Chisq)    
## mod.snow.tree.surv1               
## mod.snow.tree.surv3          1    
## mod.snow.tree.surv5     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sjp.lmer(mod.snow.tree.surv5, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.snow.tree.surv5, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
##    Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##     1308     1358     -643     1286      684 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7248 -0.6124  0.2664  0.6506  2.6048 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.09498  0.3082  
##  sp.         (Intercept) 0.56548  0.7520  
##  census      (Intercept) 1.03429  1.0170  
## Number of obs: 695, groups:  quad.unique, 60; sp., 15; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.12850    1.47414  -0.087  0.93054    
## site1                0.12602    0.11359   1.109  0.26726    
## site2               -0.02943    0.10722  -0.274  0.78371    
## treatmentS           0.48850    0.16922   2.887  0.00389 ** 
## seasonsp            -0.83093    0.84608  -0.982  0.32605    
## A.con_curt          -0.02656    0.02628  -1.011  0.31208    
## A.het_curt           0.06739    0.07128   0.946  0.34440    
## treatmentS:seasonsp -1.39301    0.22668  -6.145 7.98e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Maybe the model considering season as fixed effect makes more sense?

Survival analysis

# how about the shrubs?
#### shrub
mod.snow.shrub.surv1 <- glmer(cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt +
                                (1|quad.unique) + (1|sp.) + (1|census),
                              data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial, 
                              control=glmerControl(optimizer="optimx", 
                                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
mod.snow.shrub.surv2 <- glmer(cbind(survs, deaths) ~ site + treatment*duration + 
                                (1|quad.unique) + (1|sp.),
                              data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial, 
                              control=glmerControl(optimizer="optimx", 
                                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.snow.shrub.surv3 <- glmer(cbind(survs, deaths) ~ site + treatment*season + 
                                (1|quad.unique) + (1|sp.) + (1|census),
                              data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial, 
                              control=glmerControl(optimizer="optimx", 
                                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.snow.shrub.surv4 <- glmer(cbind(survs, deaths) ~ site + treatment + season + 
                                (1|quad.unique) + (1|sp.) + (1|census),
                              data=subset(snowdat.ag_adult, growth.form=='shrub'), family=binomial, 
                              control=glmerControl(optimizer="optimx", 
                                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.snow.shrub.surv1, mod.snow.shrub.surv2, mod.snow.shrub.surv3, mod.snow.shrub.surv4)
## Data: subset(snowdat.ag_adult, growth.form == "shrub")
## Models:
## mod.snow.shrub.surv1: cbind(survs, deaths) ~ site + treatment + A.con_curt + A.het_curt + 
## mod.snow.shrub.surv1:     (1 | quad.unique) + (1 | sp.) + (1 | census)
## mod.snow.shrub.surv2: cbind(survs, deaths) ~ site + treatment * duration + (1 | quad.unique) + 
## mod.snow.shrub.surv2:     (1 | sp.)
## mod.snow.shrub.surv4: cbind(survs, deaths) ~ site + treatment + season + (1 | quad.unique) + 
## mod.snow.shrub.surv4:     (1 | sp.) + (1 | census)
## mod.snow.shrub.surv3: cbind(survs, deaths) ~ site + treatment * season + (1 | quad.unique) + 
## mod.snow.shrub.surv3:     (1 | sp.) + (1 | census)
##                      Df    AIC    BIC  logLik deviance   Chisq Chi Df
## mod.snow.shrub.surv1  8 295.55 328.72 -139.77   279.55               
## mod.snow.shrub.surv2  8 314.37 347.54 -149.19   298.37  0.0000      0
## mod.snow.shrub.surv4  8 297.13 330.30 -140.57   281.13 17.2393      0
## mod.snow.shrub.surv3  9 298.60 335.92 -140.30   280.60  0.5299      1
##                      Pr(>Chisq)    
## mod.snow.shrub.surv1               
## mod.snow.shrub.surv2     1.0000    
## mod.snow.shrub.surv4     <2e-16 ***
## mod.snow.shrub.surv3     0.4666    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.shrub.surv.diags <- augment(mod.snow.shrub.surv3)
qplot(data=mod.snow.shrub.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qqPlot(resid(mod.snow.shrub.surv3))

summary(mod.snow.shrub.surv3, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cbind(survs, deaths) ~ site + treatment * season + (1 | quad.unique) +  
##     (1 | sp.) + (1 | census)
##    Data: subset(snowdat.ag_adult, growth.form == "shrub")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    298.6    335.9   -140.3    280.6      458 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.8448  0.0851  0.1316  0.2548  1.7041 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.4779   0.6913  
##  sp.         (Intercept) 0.6600   0.8124  
##  census      (Intercept) 1.2459   1.1162  
## Number of obs: 467, groups:  quad.unique, 52; sp., 14; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.65939    0.95875   4.860 1.17e-06 ***
## site1                0.07082    0.28253   0.251   0.8021    
## site2                0.14578    0.28483   0.512   0.6088    
## treatmentS          -1.59249    0.75305  -2.115   0.0345 *  
## seasonsp            -1.14958    1.16982  -0.983   0.3258    
## treatmentS:seasonsp -0.60741    0.81576  -0.745   0.4565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##------------------------- Density dependence-------------------------##

## For all tree species together
ggplot(data=subset(snowdat.ag_adult, growth.form=='tree'), 
       aes(x=total, y=survs/total, colour=treatment)) + 
  geom_jitter()+
  geom_smooth(aes(weight = total), method='glm', method.args=list(family="binomial")) 

# no density dependence for seedlings overall
snowdat.ag_adult$dens <- (snowdat.ag_adult$survs - mean(snowdat.ag_adult$survs))/
  (2*sd(snowdat.ag_adult$survs))

mod.snow.tree.surv_ddst1 <- glmer(cbind(survs, deaths) ~ site + treatment*season + 
                                    A.con_curt + A.het_curt + dens + 
                                    (1|quad.unique) + (1|sp.) + (1|census),
                                  data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                                  control=glmerControl(optimizer="optimx",
                                                       optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.snow.tree.surv_ddst2 <- glmer(cbind(survs, deaths) ~ treatment*duration + 
                                    A.con_curt + A.het_curt + dens + 
                                    (1|quad.unique) + (1|sp.),
                                  data=subset(snowdat.ag_adult, growth.form=='tree'), family=binomial, 
                                  control=glmerControl(optimizer="optimx",
                                                       optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.snow.tree.surv5, mod.snow.tree.surv_ddst1, mod.snow.tree.surv_ddst2)
## Data: subset(snowdat.ag_adult, growth.form == "tree")
## Models:
## mod.snow.tree.surv_ddst2: cbind(survs, deaths) ~ treatment * duration + A.con_curt + A.het_curt + 
## mod.snow.tree.surv_ddst2:     dens + (1 | quad.unique) + (1 | sp.)
## mod.snow.tree.surv5: cbind(survs, deaths) ~ site + treatment * season + A.con_curt + 
## mod.snow.tree.surv5:     A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.surv_ddst1: cbind(survs, deaths) ~ site + treatment * season + A.con_curt + 
## mod.snow.tree.surv_ddst1:     A.het_curt + dens + (1 | quad.unique) + (1 | sp.) + (1 | 
## mod.snow.tree.surv_ddst1:     census)
##                          Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.snow.tree.surv_ddst2  9 1263.2 1304.1 -622.61   1245.2              
## mod.snow.tree.surv5      11 1308.0 1358.0 -643.01   1286.0  0.000      2
## mod.snow.tree.surv_ddst1 12 1229.2 1283.7 -602.59   1205.2 80.829      1
##                          Pr(>Chisq)    
## mod.snow.tree.surv_ddst2               
## mod.snow.tree.surv5               1    
## mod.snow.tree.surv_ddst1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.treeddst.surv.diags <- augment(mod.snow.tree.surv_ddst1)
qplot(data=mod.snow.treeddst.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qqPlot(resid(mod.snow.tree.surv_ddst1))

summary(mod.snow.tree.surv_ddst1, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +  
##     A.het_curt + dens + (1 | quad.unique) + (1 | sp.) + (1 |      census)
##    Data: subset(snowdat.ag_adult, growth.form == "tree")
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1229.2   1283.7   -602.6   1205.2      683 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3732 -0.5405  0.3206  0.6915  2.4750 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.08275  0.2877  
##  sp.         (Intercept) 0.53373  0.7306  
##  census      (Intercept) 0.58770  0.7666  
## Number of obs: 695, groups:  quad.unique, 60; sp., 15; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -0.80072    1.41077  -0.568  0.57032    
## site1                0.21745    0.11167   1.947  0.05151 .  
## site2               -0.17298    0.10577  -1.635  0.10196    
## treatmentS           0.52767    0.16768   3.147  0.00165 ** 
## seasonsp            -0.50795    0.64769  -0.784  0.43290    
## A.con_curt          -0.02498    0.02609  -0.957  0.33840    
## A.het_curt           0.10022    0.07068   1.418  0.15623    
## dens                 0.88569    0.10374   8.537  < 2e-16 ***
## treatmentS:seasonsp -1.24057    0.23146  -5.360 8.33e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##------------------------- individual species -----------------------------##
## And now pulling the species apart. 
## First need to filter down to species that can be analysed - need enough
## reps (lets say 10) and variation in abundance (at least >5 
sp.snow.surv.counts <- summarise(group_by(snowdat.ag_adult, sp.), 
                                 n.quadrat=length(survs), max.dens=max(survs), min.dens=min(survs))

sp.snow.surv.counts$range <- sp.snow.surv.counts$max.dens - sp.snow.surv.counts$min.dens
sp.snow.surv.sel <- sp.snow.surv.counts$sp.[sp.snow.surv.counts$n.quadrat > 10 & sp.snow.surv.counts$range >5]
sp.snow.surv.counts[sp.snow.surv.counts$sp. %in% sp.snow.surv.sel,]
## # A tibble: 5 × 5
##      sp. n.quadrat max.dens min.dens range
##   <fctr>     <int>    <int>    <int> <int>
## 1   ABNE        17        9        0     9
## 2   ACPS        52        6        0     6
## 3   EUVE        39        6        0     6
## 4   FRMA       217       14        0    14
## 5   TIAM       189       11        0    11
ggplot(data=subset(snowdat.ag_adult, sp. %in% sp.snow.surv.sel), 
       aes(x=total, y=survs/total, colour=treatment)) + 
  geom_jitter()+
  geom_smooth(aes(weight = total), method='glm',
              method.args=list(family="binomial"), fullrange=T) +
  facet_wrap(~ sp.,  scales='free')
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

# it seems only six species have enough data, but conference intervals of EUVE are too huge,
# ABNE and ACTE are excluded becasue meeting issues when fitting models

## quick plot of only the species for which there appear to be enough data
sp.snow.surv.sel2 <- c('ACPS','FRMA','TIAM')

snow.survmods <- sapply(sp.snow.surv.sel2, function(sp){
  print(sp)
  spdat <- filter(snowdat.ag_adult, sp. == sp)
  spdat <- droplevels(spdat)
  mod <- glmer(cbind(survs, deaths) ~ site + treatment*season + dens +
                 A.con_curt + A.het_curt + (1|quad.unique) + (1|census), 
               data=spdat, family=binomial, 
               control=glmerControl(optimizer="optimx",
                                    optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
  )
}, simplify=FALSE)
## [1] "ACPS"
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
## [1] "FRMA"
## [1] "TIAM"
## diagnostic plots
indmods.s.surv.resids <- lapply(snow.survmods, augment) 

indmods.s.surv.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.s.surv.resids)), dat=indmods.s.surv.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.s.surv.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'

library(lattice)
lapply(snow.survmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACPS
## $ACPS$quad.unique

## 
## $ACPS$census

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## $FRMA$census

## 
## 
## $TIAM
## $TIAM$quad.unique

## 
## $TIAM$census

lapply(snow.survmods, function(x) summary(x, correlation=FALSE))
## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX
## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | census)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##     75.5     97.0    -26.7     53.5       41 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5213 -0.3807  0.1999  0.5332  2.0530 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  quad.unique (Intercept) 8.874e-05 0.00942 
##  census      (Intercept) 1.755e+00 1.32469 
## Number of obs: 52, groups:  quad.unique, 26; census, 6
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)           19.1402     8.9994   2.127   0.0334 *
## site2                 -1.2808     0.7986  -1.604   0.1088  
## site3                  0.9032     1.3712   0.659   0.5101  
## treatmentS            19.8869  8325.8311   0.002   0.9981  
## seasonsp              -2.4824     1.5969  -1.554   0.1201  
## dens                   0.9190     0.8284   1.109   0.2672  
## A.con_curt            -0.2111     0.6469  -0.326   0.7442  
## A.het_curt            -0.8704     0.4864  -1.790   0.0735 .
## treatmentS:seasonsp  -19.5511  8325.8311  -0.002   0.9981  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
## 
## 
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | census)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    488.5    525.7   -233.3    466.5      206 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7912 -0.6794  0.3483  0.7508  2.2526 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0878   0.2963  
##  census      (Intercept) 0.3931   0.6270  
## Number of obs: 217, groups:  quad.unique, 54; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.48326    1.86984  -0.793  0.42763    
## site2               -1.24517    0.38887  -3.202  0.00136 ** 
## site3               -0.70858    0.32464  -2.183  0.02906 *  
## treatmentS           0.51463    0.26828   1.918  0.05508 .  
## seasonsp            -0.20714    0.56601  -0.366  0.71439    
## dens                 0.90654    0.16632   5.451 5.02e-08 ***
## A.con_curt           0.05652    0.03820   1.479  0.13905    
## A.het_curt           0.14065    0.10067   1.397  0.16237    
## treatmentS:seasonsp -1.06678    0.35881  -2.973  0.00295 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## cbind(survs, deaths) ~ site + treatment * season + dens + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | census)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    408.3    444.0   -193.2    386.3      178 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.70668 -0.58270 -0.06585  0.76267  2.07260 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.07743  0.2783  
##  census      (Intercept) 0.59358  0.7704  
## Number of obs: 189, groups:  quad.unique, 59; census, 6
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.446536   2.134263   0.678    0.498    
## site2                0.130718   0.257250   0.508    0.611    
## site3               -0.079597   0.246263  -0.323    0.747    
## treatmentS           0.110885   0.231615   0.479    0.632    
## seasonsp            -1.053821   0.806889  -1.306    0.192    
## dens                 1.192432   0.206377   5.778 7.56e-09 ***
## A.con_curt          -0.140947   0.078349  -1.799    0.072 .  
## A.het_curt           0.001488   0.102478   0.015    0.988    
## treatmentS:seasonsp -0.456131   0.378106  -1.206    0.228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lapply(snow.survmods, function(x) anova(x))
## $ACPS
## Analysis of Variance Table
##                  Df Sum Sq Mean Sq F value
## site              2 2.7392  1.3696  1.3696
## treatment         1 0.7594  0.7594  0.7594
## season            1 2.8410  2.8410  2.8410
## dens              1 0.6961  0.6961  0.6961
## A.con_curt        1 0.2729  0.2729  0.2729
## A.het_curt        1 3.2023  3.2023  3.2023
## treatment:season  1 0.0000  0.0000  0.0000
## 
## $FRMA
## Analysis of Variance Table
##                  Df Sum Sq Mean Sq F value
## site              2  2.750   1.375  1.3749
## treatment         1  0.705   0.705  0.7046
## season            1  3.322   3.322  3.3221
## dens              1 32.034  32.034 32.0340
## A.con_curt        1  0.797   0.797  0.7965
## A.het_curt        1  1.471   1.471  1.4707
## treatment:season  1  8.849   8.849  8.8485
## 
## $TIAM
## Analysis of Variance Table
##                  Df Sum Sq Mean Sq F value
## site              2  0.348   0.174  0.1740
## treatment         1  0.034   0.034  0.0342
## season            1  5.277   5.277  5.2770
## dens              1 41.220  41.220 41.2204
## A.con_curt        1  3.807   3.807  3.8070
## A.het_curt        1  0.001   0.001  0.0005
## treatment:season  1  1.471   1.471  1.4712
##---------------- Survival of other tree seedlings-------------------##
# sp.sel2 <- c('ACBA', 'ACPS', 'ACTE', 'FRMA', 'PIKO', 'TIAM')
## other tree species
snowdat.ag_adult.t <- subset(snowdat.ag_adult, growth.form == 'tree')
snowdat.ag_adult.t1 <- snowdat.ag_adult.t[!(snowdat.ag_adult.t$sp. %in% sp.snow.surv.sel2),]
summary(snowdat.ag_adult.t1)
##        X          site       quadrat           sp.     treatment
##  Min.   :   5.0   1:112   Min.   :102.0   ACTE   :51   C:121    
##  1st Qu.: 144.0   2: 97   1st Qu.:203.0   ACBA   :45   S:116    
##  Median : 303.0   3: 28   Median :306.0   ACMO   :38            
##  Mean   : 349.4           Mean   :292.5   ACMA   :28            
##  3rd Qu.: 499.0           3rd Qu.:403.0   PIKO   :23            
##  Max.   :1158.0           Max.   :508.0   ABNE   :17            
##                                           (Other):35            
##  growth.form  census       deaths           total           survs       
##  shrub:  0   15sp:64   Min.   :0.0000   Min.   : 1.00   Min.   :0.0000  
##  tree :237   15fa:18   1st Qu.:0.0000   1st Qu.: 1.00   1st Qu.:0.0000  
##              16sp:12   Median :0.0000   Median : 1.00   Median :1.0000  
##              16fa:38   Mean   :0.5781   Mean   : 1.57   Mean   :0.9916  
##              17sp:33   3rd Qu.:1.0000   3rd Qu.: 2.00   3rd Qu.:1.0000  
##              17fa:72   Max.   :7.0000   Max.   :12.00   Max.   :9.0000  
##                                                                         
##   quad.unique          quad.sp        A.sum          A.con         
##  1_1_308: 18   1_1_305_ACMA:  6   Min.   :4487   Min.   :   0.000  
##  1_1_109: 12   1_1_308_ACBA:  6   1st Qu.:5335   1st Qu.:   1.712  
##  1_1_205: 10   1_1_308_ACMO:  6   Median :5657   Median :  70.242  
##  1_1_305: 10   1_1_308_ACTE:  6   Mean   :5811   Mean   : 292.147  
##  2_1_206: 10   2_1_102_ACMO:  6   3rd Qu.:6261   3rd Qu.: 315.070  
##  3_1_109: 10   2_1_104_ACTE:  6   Max.   :7822   Max.   :2670.429  
##  (Other):167   (Other)     :201                                    
##      A.het         con.dens        A.con_curt       A.het_curt   
##  Min.   :2906   Min.   : 0.000   Min.   : 0.000   Min.   :14.27  
##  1st Qu.:5085   1st Qu.: 1.000   1st Qu.: 1.196   1st Qu.:17.20  
##  Median :5450   Median : 4.000   Median : 4.126   Median :17.60  
##  Mean   :5519   Mean   : 5.738   Mean   : 4.468   Mean   :17.62  
##  3rd Qu.:6122   3rd Qu.: 8.000   3rd Qu.: 6.805   3rd Qu.:18.29  
##  Max.   :7779   Max.   :36.000   Max.   :13.874   Max.   :19.81  
##                                                                  
##    A.sum_curt    con.dens_curt      duration        season         
##  Min.   :16.49   Min.   :0.000   Min.   : 9.00   Length:237        
##  1st Qu.:17.47   1st Qu.:1.000   1st Qu.: 9.00   Class :character  
##  Median :17.82   Median :1.587   Median :24.00   Mode  :character  
##  Mean   :17.94   Mean   :1.446   Mean   :23.34                     
##  3rd Qu.:18.43   3rd Qu.:2.000   3rd Qu.:35.00                     
##  Max.   :19.85   Max.   :3.302   Max.   :35.00                     
##                                                                    
##       dens        
##  Min.   :-0.4615  
##  1st Qu.:-0.4615  
##  Median :-0.2118  
##  Mean   :-0.2139  
##  3rd Qu.:-0.2118  
##  Max.   : 1.7853  
## 
mod.snow.resttree.surv <- glmer(cbind(survs, deaths) ~ site + treatment*season + 
                                   A.con_curt + A.het_curt+
                                   (1|quad.unique) + (1|census) + (1|sp.),
                                 data=snowdat.ag_adult.t1, family=binomial, 
                                 control=glmerControl(optimizer="optimx", 
                                                      optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
## Warning in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
mod.snow.resttree.surv.diags <- augment(mod.snow.resttree.surv)
qplot(data=mod.snow.resttree.surv.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data=mod.snow.resttree.surv.diags, x=treatment, y=.wtres, geom="boxplot") 

## plot fixed and random effects
sjp.lmer(mod.snow.resttree.surv, type='fe', p.kr = FALSE)

sjp.lmer(mod.snow.resttree.surv, type='re', sort.est="sort.all", p.kr = FALSE)
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.snow.resttree.surv))

summary(mod.snow.resttree.surv, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(survs, deaths) ~ site + treatment * season + A.con_curt +  
##     A.het_curt + (1 | quad.unique) + (1 | census) + (1 | sp.)
##    Data: snowdat.ag_adult.t1
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    282.8    321.0   -130.4    260.8      226 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8257 -0.2492  0.1888  0.4405  3.4835 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.2952   0.5433  
##  sp.         (Intercept) 0.8839   0.9402  
##  census      (Intercept) 1.3861   1.1773  
## Number of obs: 237, groups:  quad.unique, 48; sp., 12; census, 6
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.9312000  4.0089492  -0.981  0.32679    
## site1                0.3240740  0.3278558   0.988  0.32292    
## site2               -0.5856986  0.3843514  -1.524  0.12754    
## treatmentS           1.7165958  0.5280091   3.251  0.00115 ** 
## seasonsp            -0.0555158  1.0901599  -0.051  0.95939    
## A.con_curt          -0.0005923  0.0979687  -0.006  0.99518    
## A.het_curt           0.2795559  0.2176052   1.285  0.19890    
## treatmentS:seasonsp -4.4869496  0.8002280  -5.607 2.06e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Parameters or bounds appear to have different scalings.
##   This can cause poor performance in optimization. 
##   It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.

Recruitment analysis

###-------------------- recruitment annlysis -----------------------------###
snow.rect.ag_adult <- read.csv("data/snow.rect.ag_adult.csv")

snow.rect.ag_adult$site <- factor(snow.rect.ag_adult$site) 
contrasts(snow.rect.ag_adult$site) <- "contr.sum"

## Trees
snow.t.rect_adult <- snow.rect.ag_adult[snow.rect.ag_adult$growth.form == 'tree',]
mod.snow.recr1 <- glmer(recruits ~ site + treatment*census + A.con_curt + A.het_curt +
                          (1|quad.unique) + (1|sp.),
                        data=snow.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
mod.snow.recr2 <- glmer(recruits ~ site + treatment + A.con_curt + A.het_curt +
                          (1|quad.unique) + (1|sp.),
                        data=snow.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))
anova(mod.snow.recr1, mod.snow.recr2)
## Data: snow.t.rect_adult
## Models:
## mod.snow.recr2: recruits ~ site + treatment + A.con_curt + A.het_curt + (1 | 
## mod.snow.recr2:     quad.unique) + (1 | sp.)
## mod.snow.recr1: recruits ~ site + treatment * census + A.con_curt + A.het_curt + 
## mod.snow.recr1:     (1 | quad.unique) + (1 | sp.)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.snow.recr2  8 1426.1 1456.2 -705.07   1410.1                         
## mod.snow.recr1 12 1348.3 1393.3 -662.16   1324.3 85.825      4  < 2.2e-16
##                   
## mod.snow.recr2    
## mod.snow.recr1 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sjp.lmer(mod.snow.recr1, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

## Random effects ofABNE, FRMA and TIAM go extramely, try to put group them against other species as the fixed effect
snow.t.rect_adult$sp.group <- ifelse(snow.t.rect_adult$sp.=='ABNE' | 
                                       snow.t.rect_adult$sp.=='FRMA' | snow.t.rect_adult$sp.=='TIAM', 1, 0)
snow.t.rect_adult$sp.group <- as.factor(snow.t.rect_adult$sp.group) 

mod.snow.recr3 <- glmer(recruits ~ site + treatment*census + 
                          A.con_curt + A.het_curt + sp.group +
                          (1|quad.unique) + (1|sp.),
                        data=snow.t.rect_adult, family=poisson, 
                        control=glmerControl(optimizer="optimx", 
                                             optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

anova(mod.snow.recr1, mod.snow.recr2, mod.snow.recr3)
## Data: snow.t.rect_adult
## Models:
## mod.snow.recr2: recruits ~ site + treatment + A.con_curt + A.het_curt + (1 | 
## mod.snow.recr2:     quad.unique) + (1 | sp.)
## mod.snow.recr1: recruits ~ site + treatment * census + A.con_curt + A.het_curt + 
## mod.snow.recr1:     (1 | quad.unique) + (1 | sp.)
## mod.snow.recr3: recruits ~ site + treatment * census + A.con_curt + A.het_curt + 
## mod.snow.recr3:     sp.group + (1 | quad.unique) + (1 | sp.)
##                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mod.snow.recr2  8 1426.1 1456.2 -705.07   1410.1                         
## mod.snow.recr1 12 1348.3 1393.3 -662.16   1324.3 85.825      4  < 2.2e-16
## mod.snow.recr3 13 1333.1 1381.9 -653.57   1307.1 17.174      1  3.412e-05
##                   
## mod.snow.recr2    
## mod.snow.recr1 ***
## mod.snow.recr3 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.recr.diags <- augment(mod.snow.recr1)
qplot(data=mod.snow.recr.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.snow.recr.diags, x=treatment, y=.wtres, geom="boxplot") 

sjp.lmer(mod.snow.recr3, type='fe')

sjp.lmer(mod.snow.recr3, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.snow.recr3))

summary(mod.snow.recr3, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +  
##     sp.group + (1 | quad.unique) + (1 | sp.)
##    Data: snow.t.rect_adult
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1333.1   1381.9   -653.6   1307.1      302 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9358 -0.7704 -0.2309  0.4938  5.0218 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.05735  0.2395  
##  sp.         (Intercept) 0.03437  0.1854  
## Number of obs: 315, groups:  quad.unique, 60; sp., 14
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.986064   0.866612   1.138  0.25519    
## site1                 -0.076718   0.070908  -1.082  0.27927    
## site2                  0.118672   0.065411   1.814  0.06964 .  
## treatmentS            -1.014233   0.329974  -3.074  0.00211 ** 
## census16sp             0.640176   0.150175   4.263 2.02e-05 ***
## census17sp             0.183148   0.164943   1.110  0.26684    
## A.con_curt             0.008625   0.014882   0.580  0.56220    
## A.het_curt            -0.062159   0.046168  -1.346  0.17818    
## sp.group1              1.018866   0.172732   5.899 3.67e-09 ***
## treatmentS:census16sp  0.911925   0.333331   2.736  0.00622 ** 
## treatmentS:census17sp  0.916593   0.343051   2.672  0.00754 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###---------------- individual species analysis ----------------------------###
## reps (lets say 10) and variation in abundance (at least >5 
sp.snowrect.counts <- summarise(group_by(snow.t.rect_adult, sp.), 
                                n.quadrat=length(recruits), max.rect=max(recruits), min.rect=min(recruits))

sp.snowrect.counts$range <- sp.snowrect.counts$max.rect - sp.snowrect.counts$min.rect
sp.snowrect.sel <- sp.snowrect.counts$sp.[sp.snowrect.counts$n.quadrat > 10 & sp.snowrect.counts$range >5]
sp.snowrect.counts[sp.snowrect.counts$sp. %in% sp.snowrect.sel,]
## # A tibble: 4 × 5
##      sp. n.quadrat max.rect min.rect range
##   <fctr>     <int>    <int>    <int> <int>
## 1   ABNE        16       12        1    11
## 2   ACPS        30        9        1     8
## 3   FRMA        99       19        1    18
## 4   TIAM       104       12        1    11
sp.snowrect.sel2 <- c('ACPS','FRMA','TIAM')

snowrectmods <- sapply(sp.snowrect.sel2, function(sp){
  print(sp)
  spdat <- filter(snow.t.rect_adult, sp. == sp)
  spdat <- droplevels(spdat)
  mod <- glmer(recruits ~ site + treatment*census +
                 A.con_curt + A.het_curt + (1|quad.unique), 
               data=spdat, family=poisson, 
               control=glmerControl(optimizer="optimx",
                                    optCtrl=list(method=c('bobyqa', 'Nelder-Mead')))
  )
}, simplify=FALSE)
## [1] "ACPS"
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
## [1] "FRMA"
## [1] "TIAM"
## diagnostic plots
indmods.snowrect.resids <- lapply(snowrectmods, augment) 

indmods.snowrect.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.snowrect.resids)), dat=indmods.snowrect.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.snowrect.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : span too small. fewer data values than degrees of freedom.
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -0.0033838
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.059818
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.43823
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small.
## fewer data values than degrees of freedom.
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -0.0033838
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.059818
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 0.43823

lapply(snowrectmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACPS
## $ACPS$quad.unique

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## 
## $TIAM
## $TIAM$quad.unique

lapply(snowrectmods, function(x) summary(x, correlation=FALSE))
## $ACPS
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +  
##     (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    112.2    124.8    -47.1     94.2       21 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3579 -0.5173 -0.1187  0.3802  2.6937 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  quad.unique (Intercept) 0.0003189 0.01786 
## Number of obs: 30, groups:  quad.unique, 24
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -3.5180     3.7233  -0.945    0.345
## site2         0.3474     0.3244   1.071    0.284
## site3        -0.3199     0.5520  -0.580    0.562
## treatmentS   -0.1368     0.6354  -0.215    0.830
## census16sp    0.4643     1.1322   0.410    0.682
## census17sp    1.0981     1.1792   0.931    0.352
## A.con_curt    0.2011     0.2257   0.891    0.373
## A.het_curt    0.1149     0.1944   0.591    0.554
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
## 
## 
## $FRMA
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +  
##     (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    438.8    467.3   -208.4    416.8       88 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7627 -0.8554 -0.1045  0.6032  3.0342 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.02111  0.1453  
## Number of obs: 99, groups:  quad.unique, 54
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            2.08221    1.03230   2.017 0.043690 *  
## site2                  0.99959    0.18393   5.435 5.49e-08 ***
## site3                  0.45709    0.16984   2.691 0.007118 ** 
## treatmentS            -1.13180    0.38565  -2.935 0.003338 ** 
## census16sp             0.62764    0.16391   3.829 0.000129 ***
## census17sp            -0.73662    0.33055  -2.228 0.025848 *  
## A.con_curt            -0.06822    0.01976  -3.452 0.000556 ***
## A.het_curt            -0.06138    0.05608  -1.094 0.273769    
## treatmentS:census16sp  1.09505    0.40032   2.735 0.006229 ** 
## treatmentS:census17sp  1.28666    0.51631   2.492 0.012701 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $TIAM
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +  
##     (1 | quad.unique)
##    Data: spdat
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    455.3    484.4   -216.7    433.3       93 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2675 -0.5815 -0.1436  0.3405  2.2702 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.1564   0.3954  
## Number of obs: 104, groups:  quad.unique, 58
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)   
## (Intercept)           -0.920613   1.667807  -0.552  0.58096   
## site2                 -0.103781   0.194799  -0.533  0.59420   
## site3                 -0.017650   0.194939  -0.090  0.92786   
## treatmentS             0.041453   0.905515   0.046  0.96349   
## census16sp             1.509932   0.529904   2.849  0.00438 **
## census17sp             0.852936   0.538295   1.584  0.11308   
## A.con_curt             0.087065   0.061449   1.417  0.15652   
## A.het_curt            -0.002372   0.078127  -0.030  0.97578   
## treatmentS:census16sp -0.175931   0.910975  -0.193  0.84686   
## treatmentS:census17sp  0.267613   0.919976   0.291  0.77113   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lapply(snowrectmods, function(x) anova(x))
## $ACPS
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 1.47123 0.73561  0.7356
## treatment   1 1.97066 1.97066  1.9707
## census      2 2.85826 1.42913  1.4291
## A.con_curt  1 0.89657 0.89657  0.8966
## A.het_curt  1 0.34909 0.34909  0.3491
## 
## $FRMA
## Analysis of Variance Table
##                  Df Sum Sq Mean Sq F value
## site              2 13.021   6.510  6.5105
## treatment         1  0.741   0.741  0.7409
## census            2 72.516  36.258 36.2582
## A.con_curt        1 10.690  10.690 10.6902
## A.het_curt        1  1.313   1.313  1.3128
## treatment:census  2  8.029   4.014  4.0144
## 
## $TIAM
## Analysis of Variance Table
##                  Df  Sum Sq Mean Sq F value
## site              2  1.3924  0.6962  0.6962
## treatment         1  0.0314  0.0314  0.0314
## census            2 22.5802 11.2901 11.2901
## A.con_curt        1  2.3351  2.3351  2.3351
## A.het_curt        1  0.0000  0.0000  0.0000
## treatment:census  2  3.9523  1.9762  1.9762
###----------------- Pooling other species -------------------------###

snow.t.rect.other_adult <- snow.t.rect_adult[!(snow.t.rect_adult$sp. %in% sp.snowrect.sel2),]
summary(snow.t.rect.other_adult)
##        X          site      quadrat           sp.     treatment
##  Min.   :  3.00   1:32   Min.   :102.0   ABNE   :16   C:30     
##  1st Qu.: 81.25   2:37   1st Qu.:203.0   ACBA   :15   S:52     
##  Median :142.00   3:13   Median :305.0   ACMO   :12            
##  Mean   :153.76          Mean   :290.1   QUMO   :10            
##  3rd Qu.:225.75          3rd Qu.:404.0   ACTE   : 8            
##  Max.   :343.00          Max.   :508.0   PIKO   : 7            
##                                          (Other):14            
##  growth.form  census      recruits       quad.unique         quad.sp  
##  shrub: 0    15sp: 2   Min.   : 1.000   1_1_205: 5   1_1_205_ACBA: 2  
##  tree :82    16sp:30   1st Qu.: 1.000   1_1_508: 5   1_1_409_ACBA: 2  
##              17sp:50   Median : 1.000   1_1_109: 4   1_1_508_ACBA: 2  
##                        Mean   : 1.829   1_1_502: 4   2_1_203_ABNE: 2  
##                        3rd Qu.: 2.000   2_1_203: 4   2_1_501_ACTE: 2  
##                        Max.   :12.000   3_1_109: 4   3_1_109_ACBA: 2  
##                                         (Other):56   (Other)     :70  
##      A.sum          A.con             A.het         con.dens     
##  Min.   :4487   Min.   :   0.00   Min.   :2906   Min.   : 0.000  
##  1st Qu.:5154   1st Qu.:  15.69   1st Qu.:4530   1st Qu.: 1.000  
##  Median :5627   Median : 191.49   Median :5267   Median : 5.500  
##  Mean   :5699   Mean   : 427.85   Mean   :5271   Mean   : 8.451  
##  3rd Qu.:6261   3rd Qu.: 353.46   3rd Qu.:5955   3rd Qu.: 9.750  
##  Max.   :7822   Max.   :2670.43   Max.   :7465   Max.   :36.000  
##                                                                  
##    A.con_curt       A.het_curt      A.sum_curt    con.dens_curt   sp.group
##  Min.   : 0.000   Min.   :14.27   Min.   :16.49   Min.   :0.000   0:66    
##  1st Qu.: 2.489   1st Qu.:16.55   1st Qu.:17.27   1st Qu.:1.000   1:16    
##  Median : 5.764   Median :17.40   Median :17.79   Median :1.764           
##  Mean   : 5.538   Mean   :17.34   Mean   :17.83   Mean   :1.672           
##  3rd Qu.: 7.070   3rd Qu.:18.13   3rd Qu.:18.43   3rd Qu.:2.136           
##  Max.   :13.874   Max.   :19.54   Max.   :19.85   Max.   :3.302           
## 
mod.snow.recr.other <- glmer(recruits ~ site + treatment*census + A.con_curt + A.het_curt +
                                (1|quad.unique) + (1|sp.),
                              data=snow.t.rect.other_adult, family=poisson, 
                              control=glmerControl(optimizer="optimx", 
                                                   optCtrl=list(method=c('bobyqa', 'Nelder-Mead'))))

mod.snow.recr.other.diags <- augment(mod.snow.recr.other)
qplot(data=mod.snow.recr.other.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.snow.recr.other.diags, x=treatment, y=.wtres, geom="boxplot") 

qqPlot(resid(mod.snow.recr.other))

sjp.lmer(mod.snow.recr.other, type='fe')

sjp.lmer(mod.snow.recr.other, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.snow.recr.other, correlation=FALSE)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: recruits ~ site + treatment * census + A.con_curt + A.het_curt +  
##     (1 | quad.unique) + (1 | sp.)
##    Data: snow.t.rect.other_adult
## Control: 
## glmerControl(optimizer = "optimx", optCtrl = list(method = c("bobyqa",  
##     "Nelder-Mead")))
## 
##      AIC      BIC   logLik deviance df.resid 
##    268.3    297.2   -122.2    244.3       70 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9176 -0.3487 -0.1552  0.1270  3.4479 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.005358 0.0732  
##  sp.         (Intercept) 0.157223 0.3965  
## Number of obs: 82, groups:  quad.unique, 42; sp., 11
## 
## Fixed effects:
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)            0.93467    1.90000   0.492    0.623
## site1                  0.08930    0.16227   0.550    0.582
## site2                  0.04132    0.20385   0.203    0.839
## treatmentS             0.11545    1.41667   0.082    0.935
## census16sp             1.08750    1.06984   1.016    0.309
## census17sp             1.62865    1.01503   1.604    0.109
## A.con_curt             0.00336    0.04263   0.079    0.937
## A.het_curt            -0.10866    0.08974  -1.211    0.226
## treatmentS:census16sp -0.16913    1.47120  -0.115    0.908
## treatmentS:census17sp -0.44983    1.43142  -0.314    0.753

Growth analysis

###----------------------- Growth models ------------------------------###
snowdat.g_adult <- read.csv("data/snowdat.g_adult.csv")
snowdat.g_adult$site <- factor(snowdat.g_adult$site)
contrasts(snowdat.g_adult$site) <- "contr.sum"

#### duration
## duration is presented as months from the beginning of the experiment
snowdat.g_adult$duration <- ifelse(snowdat.g_adult$census=="15sp", 9,
                                   ifelse(snowdat.g_adult$census=="15fa", 12,
                                          ifelse(snowdat.g_adult$census=="16sp", 21,
                                                 ifelse(snowdat.g_adult$census=="16fa", 24,
                                                        ifelse(snowdat.g_adult$census=="17sp", 32,
                                                               ifelse(snowdat.g_adult$census=="17fa", 35, NA))))))
## season
snowdat.g_adult$season <- ifelse(snowdat.g_adult$census=='15fa'| 
                                   snowdat.g_adult$census=='16fa'|
                                   snowdat.g_adult$census=='17fa', 'fa', 'sp')

summary(snowdat.g_adult)
##        X             sp.      site       quadrat           tag       
##  Min.   :  81   FRMA   :793   1:695   Min.   :102.0   Min.   : 1001  
##  1st Qu.:5138   TIAM   :379   2:926   1st Qu.:203.0   1st Qu.: 3008  
##  Median :6489   PHSC   :179   3:462   Median :306.0   Median : 5031  
##  Mean   :6114   EUMA   : 96           Mean   :295.9   Mean   : 5513  
##  3rd Qu.:8017   ACPS   : 75           3rd Qu.:404.0   3rd Qu.: 8020  
##  Max.   :9210   EUVE   : 69           Max.   :508.0   Max.   :10051  
##                 (Other):492                                          
##  treatment growth.form   census       grow.snd         quad.unique  
##  C:1211    shrub: 612   15fa:142   Min.   :-2.48608   2_1_210:  88  
##  S: 872    tree :1471   15sp:120   1st Qu.:-0.42166   2_1_403:  83  
##                         16fa:589   Median : 0.08960   2_1_501:  76  
##                         16sp:117   Mean   : 0.03001   1_1_409:  69  
##                         17fa:693   3rd Qu.: 0.44315   2_1_104:  68  
##                         17sp:422   Max.   : 2.53909   2_1_401:  67  
##                                                       (Other):1632  
##          quad.sp         A.sum          A.con            A.het     
##  2_1_210_FRMA:  50   Min.   :4487   Min.   :   0.0   Min.   :2827  
##  2_1_102_FRMA:  38   1st Qu.:5335   1st Qu.:   0.0   1st Qu.:4573  
##  2_1_403_FRMA:  35   Median :5708   Median : 258.6   Median :5373  
##  2_1_507_FRMA:  34   Mean   :5899   Mean   : 577.9   Mean   :5321  
##  2_1_308_FRMA:  33   3rd Qu.:6533   3rd Qu.: 981.2   3rd Qu.:5974  
##  2_1_401_FRMA:  28   Max.   :9223   Max.   :4003.9   Max.   :9223  
##  (Other)     :1865                                                 
##     con.dens        A.con_curt       A.het_curt      A.sum_curt   
##  Min.   : 0.000   Min.   : 0.000   Min.   :14.14   Min.   :16.49  
##  1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.:16.60   1st Qu.:17.47  
##  Median : 2.000   Median : 6.371   Median :17.51   Median :17.87  
##  Mean   : 5.039   Mean   : 5.469   Mean   :17.38   Mean   :18.02  
##  3rd Qu.: 8.000   3rd Qu.: 9.937   3rd Qu.:18.14   3rd Qu.:18.69  
##  Max.   :41.000   Max.   :15.879   Max.   :20.97   Max.   :20.97  
##                                                                   
##  con.dens_curt      duration        season         
##  Min.   :0.000   Min.   : 9.00   Length:2083       
##  1st Qu.:0.000   1st Qu.:24.00   Class :character  
##  Median :1.260   Median :32.00   Mode  :character  
##  Mean   :1.114   Mean   :27.43                     
##  3rd Qu.:2.000   3rd Qu.:35.00                     
##  Max.   :3.448   Max.   :35.00                     
## 
## tree
mod.snow.tree.grow1 <- lmer(grow.snd ~ site + treatment + A.con_curt + A.het_curt +
                              (1|quad.unique) + (1|census) + (1|sp.), 
                            data=subset(snowdat.g_adult, growth.form == 'tree'))

mod.snow.tree.grow2 <- lmer(grow.snd ~ site + treatment*duration+
                              (1|quad.unique) + (1|sp.), 
                            data=subset(snowdat.g_adult, growth.form == 'tree'))

mod.snow.tree.grow3 <- lmer(grow.snd ~ site + treatment*season+
                              (1|quad.unique) + (1|sp.) + (1|census), 
                            data=subset(snowdat.g_adult, growth.form == 'tree'))

anova(mod.snow.tree.grow1, mod.snow.tree.grow2, mod.snow.tree.grow3)
## refitting model(s) with ML (instead of REML)
## Data: subset(snowdat.g_adult, growth.form == "tree")
## Models:
## mod.snow.tree.grow2: grow.snd ~ site + treatment * duration + (1 | quad.unique) + 
## mod.snow.tree.grow2:     (1 | sp.)
## mod.snow.tree.grow1: grow.snd ~ site + treatment + A.con_curt + A.het_curt + (1 | 
## mod.snow.tree.grow1:     quad.unique) + (1 | census) + (1 | sp.)
## mod.snow.tree.grow3: grow.snd ~ site + treatment * season + (1 | quad.unique) + (1 | 
## mod.snow.tree.grow3:     sp.) + (1 | census)
##                     Df    AIC    BIC  logLik deviance    Chisq Chi Df
## mod.snow.tree.grow2  9 3624.0 3671.7 -1803.0   3606.0                
## mod.snow.tree.grow1 10 3433.3 3486.3 -1706.7   3413.3 192.7130      1
## mod.snow.tree.grow3 10 3432.7 3485.6 -1706.3   3412.7   0.6714      0
##                     Pr(>Chisq)    
## mod.snow.tree.grow2               
## mod.snow.tree.grow1  < 2.2e-16 ***
## mod.snow.tree.grow3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.tree.grow.diags <- augment(mod.snow.tree.grow1)
qplot(data=mod.snow.tree.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'gam'

qplot(data = mod.snow.tree.grow.diags, x=treatment, y=.wtres, geom="boxplot") 

qqPlot(resid(mod.snow.tree.grow1))

summary(mod.snow.tree.grow1, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.con_curt + A.het_curt + (1 |  
##     quad.unique) + (1 | census) + (1 | sp.)
##    Data: subset(snowdat.g_adult, growth.form == "tree")
## 
## REML criterion at convergence: 3444.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8592 -0.5555 -0.0079  0.5554  3.3538 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.000000 0.0000  
##  sp.         (Intercept) 0.004665 0.0683  
##  census      (Intercept) 0.142553 0.3776  
##  Residual                0.588307 0.7670  
## Number of obs: 1471, groups:  quad.unique, 59; sp., 13; census, 6
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) -0.6298832  0.4244555  -1.484
## site1        0.1347980  0.0340059   3.964
## site2       -0.0168361  0.0313303  -0.537
## treatmentS   0.0436222  0.0417431   1.045
## A.con_curt   0.0005533  0.0072676   0.076
## A.het_curt   0.0333985  0.0211320   1.580
# although treating site as a fixed effect make the model performing better, results do not change 
# duration has no effect on seedling growth, there is no need to consider it in the following analyses.


## shrub
mod.snow.shrub.grow1 <- lmer(grow.snd ~ site + treatment + A.sum_curt+
                               (1|quad.unique) + (1|census) + (1|sp.), 
                             data=subset(snowdat.g_adult, growth.form == 'shrub'))
mod.snow.shrub.grow2 <- lmer(grow.snd ~ site + treatment + 
                               (1|quad.unique) + (1|census) + (1|sp.), 
                             data=subset(snowdat.g_adult, growth.form == 'shrub'))
mod.snow.shrub.grow3 <- lmer(grow.snd ~ treatment*season + 
                               (1|quad.unique) + (1|census) + (1|sp.), 
                             data=subset(snowdat.g_adult, growth.form == 'shrub'))

anova(mod.snow.shrub.grow1, mod.snow.shrub.grow2, mod.snow.shrub.grow3)
## refitting model(s) with ML (instead of REML)
## Data: subset(snowdat.g_adult, growth.form == "shrub")
## Models:
## mod.snow.shrub.grow2: grow.snd ~ site + treatment + (1 | quad.unique) + (1 | census) + 
## mod.snow.shrub.grow2:     (1 | sp.)
## mod.snow.shrub.grow3: grow.snd ~ treatment * season + (1 | quad.unique) + (1 | census) + 
## mod.snow.shrub.grow3:     (1 | sp.)
## mod.snow.shrub.grow1: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) + 
## mod.snow.shrub.grow1:     (1 | census) + (1 | sp.)
##                      Df    AIC    BIC  logLik deviance  Chisq Chi Df
## mod.snow.shrub.grow2  8 1240.3 1275.6 -612.14   1224.3              
## mod.snow.shrub.grow3  8 1239.5 1274.8 -611.76   1223.5 0.7582      0
## mod.snow.shrub.grow1  9 1237.5 1277.2 -609.74   1219.5 4.0374      1
##                      Pr(>Chisq)    
## mod.snow.shrub.grow2               
## mod.snow.shrub.grow3     <2e-16 ***
## mod.snow.shrub.grow1     0.0445 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.snow.shrub.grow.diags <- augment(mod.snow.shrub.grow1)
qplot(data=mod.snow.shrub.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.snow.shrub.grow.diags, x=treatment, y=.wtres, geom="boxplot") 

qqPlot(resid(mod.snow.shrub.grow1))

summary(mod.snow.shrub.grow1, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census) + (1 | sp.)
##    Data: subset(snowdat.g_adult, growth.form == "shrub")
## 
## REML criterion at convergence: 1240.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7576 -0.4549  0.0047  0.4997  3.1070 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.00000  0.0000  
##  sp.         (Intercept) 0.00000  0.0000  
##  census      (Intercept) 0.05456  0.2336  
##  Residual                0.42203  0.6496  
## Number of obs: 612, groups:  quad.unique, 51; sp., 13; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  1.36225    0.60169   2.264
## site1       -0.01962    0.04158  -0.472
## site2        0.02836    0.03670   0.773
## treatmentS  -0.02627    0.05697  -0.461
## A.sum_curt  -0.07138    0.03263  -2.188
##------------------------- individual species -------------------------------##
## And now pulling the species apart. 
## First need to filter down to species that can be analysed - need enough
## reps (lets say 50) and variation in abundance (at least >10 
sp.snowgrow.counts <- summarise(group_by(snowdat.g_adult, sp.), 
                                n.indiv=length(grow.snd), n.quadrat=length(unique(grow.snd)))

sp.snowgrow.sel <- sp.snowgrow.counts$sp.[sp.snowgrow.counts$n.indiv > 50 & sp.snowgrow.counts$n.quadrat > 10]
sp.snowgrow.counts[sp.snowgrow.counts$sp. %in% sp.snowgrow.sel,]
## # A tibble: 8 × 3
##      sp. n.indiv n.quadrat
##   <fctr>   <int>     <int>
## 1   ACPS      75        14
## 2   DEGL      68        24
## 3   EUAL      65        19
## 4   EUMA      96        26
## 5   EUVE      69        14
## 6   FRMA     793        24
## 7   PHSC     179        43
## 8   TIAM     379        13
sp.snowgrow.sel2 <- c('ACPS','DEGL','EUAL','EUMA', 'EUVE','FRMA','PHSC','TIAM')

snowgrowmods <- sapply(sp.snowgrow.sel2, function(sp){
  print(sp)
  spdat <- filter(snowdat.g_adult, sp. == sp)
  spdat <- droplevels(spdat)
  mod <- lmer(grow.snd ~ site+ treatment + A.sum_curt + (1|quad.unique) + (1|census), 
              data=spdat)
}, simplify=FALSE)
## [1] "ACPS"
## [1] "DEGL"
## [1] "EUAL"
## [1] "EUMA"
## [1] "EUVE"
## [1] "FRMA"
## [1] "PHSC"
## [1] "TIAM"
## diagnostic plots
indmods.snowgrow.resids <- lapply(snowgrowmods, augment) 

indmods.snowgrow.resids <- do.call('rbind', mapply(function(spnm, dat) {
  dat$species <- spnm
  return(dat)
}, spnm=as.list(names(indmods.snowgrow.resids)), dat=indmods.snowgrow.resids, 
SIMPLIFY=FALSE))

ggplot(data=indmods.snowgrow.resids, aes(x=.fitted, y=.wtres, colour=treatment)) +
  geom_jitter() + geom_smooth() + facet_wrap(~species, scales='free') +
  geom_hline(yintercept=0, linetype='dashed')
## `geom_smooth()` using method = 'loess'

library(lattice)
lapply(snowgrowmods, function(mod){
  dotplot(ranef(mod, condVar=TRUE))})
## $ACPS
## $ACPS$quad.unique

## 
## $ACPS$census

## 
## 
## $DEGL
## $DEGL$quad.unique

## 
## $DEGL$census

## 
## 
## $EUAL
## $EUAL$quad.unique

## 
## $EUAL$census

## 
## 
## $EUMA
## $EUMA$quad.unique

## 
## $EUMA$census

## 
## 
## $EUVE
## $EUVE$quad.unique

## 
## $EUVE$census

## 
## 
## $FRMA
## $FRMA$quad.unique

## 
## $FRMA$census

## 
## 
## $PHSC
## $PHSC$quad.unique

## 
## $PHSC$census

## 
## 
## $TIAM
## $TIAM$quad.unique

## 
## $TIAM$census

lapply(snowgrowmods, function(x) summary(x, correlation=FALSE))
## $ACPS
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 157
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3769 -0.5956 -0.1252  0.7792  2.0944 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 1.518e-16 1.232e-08
##  census      (Intercept) 2.306e-01 4.802e-01
##  Residual                4.086e-01 6.392e-01
## Number of obs: 75, groups:  quad.unique, 23; census, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.92470    2.19067   0.422
## site2       -0.26127    0.20462  -1.277
## site3       -0.46614    0.20549  -2.268
## treatmentS   0.10123    0.23310   0.434
## A.sum_curt  -0.02566    0.12401  -0.207
## 
## $DEGL
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 141.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8994 -0.3190 -0.0513  0.4379  2.2436 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  census      (Intercept) 0.1721   0.4149  
##  Residual                0.3839   0.6196  
## Number of obs: 68, groups:  quad.unique, 10; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  3.86927    2.92275   1.324
## site2        0.32105    0.28425   1.129
## site3       -0.13271    0.39301  -0.338
## treatmentS   0.02295    0.17033   0.135
## A.sum_curt  -0.21884    0.16525  -1.324
## 
## $EUAL
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 153.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.63460 -0.38641 -0.04121  0.64050  2.06026 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 1.717e-15 4.144e-08
##  census      (Intercept) 1.073e-01 3.276e-01
##  Residual                5.709e-01 7.556e-01
## Number of obs: 65, groups:  quad.unique, 9; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  3.01245    2.67990   1.124
## site2       -0.01935    0.27799  -0.070
## site3       -0.37808    0.38383  -0.985
## treatmentS  -0.03388    0.79567  -0.043
## A.sum_curt  -0.15730    0.15082  -1.043
## 
## $EUMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 190.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22987 -0.45917 -0.03787  0.54778  2.72495 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 6.175e-18 2.485e-09
##  census      (Intercept) 7.848e-02 2.801e-01
##  Residual                3.712e-01 6.093e-01
## Number of obs: 96, groups:  quad.unique, 17; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  4.20258    1.67419   2.510
## site2        0.24947    0.17292   1.443
## site3        0.25304    0.28763   0.880
## treatmentS  -0.11429    0.18321  -0.624
## A.sum_curt  -0.23330    0.09523  -2.450
## 
## $EUVE
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 64.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.52899 -0.72941  0.00044  0.55533  2.95673 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  census      (Intercept) 0.3640   0.6033  
##  Residual                0.1102   0.3320  
## Number of obs: 69, groups:  quad.unique, 25; census, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -1.86794    1.52112  -1.228
## site2       -0.25097    0.13524  -1.856
## site3       -0.39067    0.13243  -2.950
## treatmentS   0.11792    0.13009   0.906
## A.sum_curt   0.09580    0.08444   1.135
## 
## $FRMA
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 1743
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9285 -0.6169  0.0257  0.5867  2.9782 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.0000   0.0000  
##  census      (Intercept) 0.2182   0.4671  
##  Residual                0.5048   0.7105  
## Number of obs: 793, groups:  quad.unique, 52; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.20195    0.51213   0.394
## site2       -0.15881    0.06735  -2.358
## site3       -0.14253    0.07743  -1.841
## treatmentS   0.06777    0.05315   1.275
## A.sum_curt  -0.01362    0.02699  -0.505
## 
## $PHSC
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 351.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6747 -0.3637 -0.0267  0.4419  3.1665 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  quad.unique (Intercept) 0.0009391 0.03064 
##  census      (Intercept) 0.0295244 0.17183 
##  Residual                0.3786474 0.61534 
## Number of obs: 179, groups:  quad.unique, 20; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.73684    0.90776   0.812
## site2        0.08695    0.10975   0.792
## site3        0.10830    0.14705   0.736
## treatmentS   0.04211    0.10331   0.408
## A.sum_curt  -0.04240    0.05140  -0.825
## 
## $TIAM
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census)
##    Data: spdat
## 
## REML criterion at convergence: 1000.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9642 -0.6932  0.1655  0.5997  2.6766 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev. 
##  quad.unique (Intercept) 8.748e-16 2.958e-08
##  census      (Intercept) 1.436e-02 1.198e-01
##  Residual                7.914e-01 8.896e-01
## Number of obs: 379, groups:  quad.unique, 55; census, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.52974    0.98108  -0.540
## site2       -0.15713    0.11993  -1.310
## site3       -0.37056    0.11569  -3.203
## treatmentS  -0.01124    0.09531  -0.118
## A.sum_curt   0.03878    0.05464   0.710
lapply(snowgrowmods, function(x) anova(x))
## $ACPS
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 3.6090 1.80449  4.4167
## treatment   1 0.0614 0.06141  0.1503
## A.sum_curt  1 0.0175 0.01750  0.0428
## 
## $DEGL
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 1.22278 0.61139  1.5925
## treatment   1 0.17727 0.17727  0.4617
## A.sum_curt  1 0.67330 0.67330  1.7538
## 
## $EUAL
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 0.56452 0.28226  0.4944
## treatment   1 0.03197 0.03197  0.0560
## A.sum_curt  1 0.62099 0.62099  1.0878
## 
## $EUMA
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 0.06733 0.03367  0.0907
## treatment   1 0.60827 0.60827  1.6385
## A.sum_curt  1 2.22815 2.22815  6.0019
## 
## $EUVE
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## site        2 0.86165 0.43083  3.9078
## treatment   1 0.04710 0.04710  0.4272
## A.sum_curt  1 0.14190 0.14190  1.2871
## 
## $FRMA
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 3.9208 1.96039  3.8832
## treatment   1 0.7907 0.79071  1.5663
## A.sum_curt  1 0.1286 0.12858  0.2547
## 
## $PHSC
## Analysis of Variance Table
##            Df  Sum Sq  Mean Sq F value
## site        2 0.17711 0.088552  0.2339
## treatment   1 0.12116 0.121164  0.3200
## A.sum_curt  1 0.25772 0.257717  0.6806
## 
## $TIAM
## Analysis of Variance Table
##            Df Sum Sq Mean Sq F value
## site        2 7.7853  3.8927  4.9186
## treatment   1 0.0702  0.0702  0.0888
## A.sum_curt  1 0.3985  0.3985  0.5036
##------------------------ Other species ------------------------##
snowdat.g.other_adult <- snowdat.g_adult[!(snowdat.g_adult$sp. %in% sp.snowgrow.sel2),]
mod.snow.other.grow <- lmer(grow.snd ~ site + treatment+  A.sum_curt +
                               (1|quad.unique) + (1|census)+ (1|sp.), data=snowdat.g.other_adult)

mod.snow.other.grow.diags <- augment(mod.snow.other.grow)
qplot(data=mod.snow.other.grow.diags, x=.fitted, y=.resid, geom=c('point', 'smooth')) +
  geom_hline(yintercept=0, linetype='dotted')
## `geom_smooth()` using method = 'loess'

qplot(data = mod.snow.other.grow.diags, x=treatment, y=.wtres, geom="boxplot") 

qqPlot(resid(mod.snow.other.grow))

summary(mod.snow.other.grow, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: grow.snd ~ site + treatment + A.sum_curt + (1 | quad.unique) +  
##     (1 | census) + (1 | sp.)
##    Data: snowdat.g.other_adult
## 
## REML criterion at convergence: 797.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4727 -0.5532  0.0165  0.5332  3.0289 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.00232  0.04817 
##  sp.         (Intercept) 0.00000  0.00000 
##  census      (Intercept) 0.06168  0.24836 
##  Residual                0.50184  0.70840 
## Number of obs: 359, groups:  quad.unique, 45; sp., 18; census, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.40890    0.97006   0.422
## site1        0.01117    0.05755   0.194
## site2       -0.04173    0.05384  -0.775
## treatmentS  -0.07238    0.08630  -0.839
## A.sum_curt  -0.01542    0.05285  -0.292

Snow removal did not affect the seedling growth.

Diversity analysis

library(vegan)
## Loading required package: permute
## This is vegan 2.4-4
divers_woodsnow <- summarise(
group_by(snowdat.ag_adult, site, quadrat, quad.unique, treatment, census), 
shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))

adult.ngbr.inf <- read.csv("data/adult.ngbr.inf.csv")
divers_woodsnow$A.sum <- adult.ngbr.inf$A.sum[match(divers_woodsnow$quad.unique, adult.ngbr.inf$quad.unique)]
summary(divers_woodsnow)
##  site       quadrat       quad.unique  treatment  census  
##  1:111   Min.   :102.0   1_1_102:  6   C:172     15sp:53  
##  2:114   1st Qu.:203.0   1_1_104:  6   S:147     15fa:46  
##  3: 94   Median :306.0   1_1_203:  6             16sp:43  
##          Mean   :299.5   1_1_205:  6             16fa:59  
##          3rd Qu.:404.0   1_1_206:  6             17sp:58  
##          Max.   :508.0   1_1_210:  6             17fa:60  
##                          (Other):283                      
##     shannon          simpson         species          A.sum     
##  Min.   :0.0000   Min.   :1.000   Min.   :0.000   Min.   :4487  
##  1st Qu.:0.0000   1st Qu.:1.684   1st Qu.:1.000   1st Qu.:5321  
##  Median :0.8018   Median :2.286   Median :3.000   Median :5662  
##  Mean   :0.8094   Mean   :  Inf   Mean   :3.028   Mean   :5900  
##  3rd Qu.:1.2730   3rd Qu.:3.341   3rd Qu.:4.000   3rd Qu.:6504  
##  Max.   :2.1458   Max.   :  Inf   Max.   :9.000   Max.   :9223  
## 
divers_woodsnow$A.sum_curt <- divers_woodsnow$A.sum^(1/3)


# season
divers_woodsnow$season <- ifelse(divers_woodsnow$census=='15fa'| 
                                   divers_woodsnow$census=='16fa'|
                                   divers_woodsnow$census=='17fa', 'fa', 'sp')
summary(divers_woodsnow)
##  site       quadrat       quad.unique  treatment  census  
##  1:111   Min.   :102.0   1_1_102:  6   C:172     15sp:53  
##  2:114   1st Qu.:203.0   1_1_104:  6   S:147     15fa:46  
##  3: 94   Median :306.0   1_1_203:  6             16sp:43  
##          Mean   :299.5   1_1_205:  6             16fa:59  
##          3rd Qu.:404.0   1_1_206:  6             17sp:58  
##          Max.   :508.0   1_1_210:  6             17fa:60  
##                          (Other):283                      
##     shannon          simpson         species          A.sum     
##  Min.   :0.0000   Min.   :1.000   Min.   :0.000   Min.   :4487  
##  1st Qu.:0.0000   1st Qu.:1.684   1st Qu.:1.000   1st Qu.:5321  
##  Median :0.8018   Median :2.286   Median :3.000   Median :5662  
##  Mean   :0.8094   Mean   :  Inf   Mean   :3.028   Mean   :5900  
##  3rd Qu.:1.2730   3rd Qu.:3.341   3rd Qu.:4.000   3rd Qu.:6504  
##  Max.   :2.1458   Max.   :  Inf   Max.   :9.000   Max.   :9223  
##                                                                 
##    A.sum_curt       season         
##  Min.   :16.49   Length:319        
##  1st Qu.:17.46   Class :character  
##  Median :17.82   Mode  :character  
##  Mean   :18.02                     
##  3rd Qu.:18.67                     
##  Max.   :20.97                     
## 
ggplot(divers_woodsnow, aes(x=treatment, y=shannon)) +  geom_boxplot()

summarise(group_by(divers_woodsnow, treatment, census), meanH=mean(shannon), 
seH=sd(shannon)/sqrt(length(shannon))) %>%
ggplot(aes(x=treatment, y=meanH, ymin=meanH-seH, ymax=meanH+seH, colour=census)) +
geom_point(position=position_dodge(width = .5))+
geom_errorbar(position=position_dodge(width = .5), width = 0.2) + 
coord_trans(y="exp") + 
labs(x="Snow removal", y="Effective number of species") 

###-------------------- Overall ---------------------------###
# shannon diversity
mod.wsnow.shan1 <- lmer(shannon ~ site + treatment + (1|quad.unique) + (1|census), 
                           data=divers_woodsnow)
mod.wsnow.shan2 <- lmer(shannon ~ site +treatment +A.sum_curt +(1|quad.unique)+ (1|census), 
                           data=divers_woodsnow)
mod.wsnow.shan3 <- lmer(shannon ~ treatment*season + (1|quad.unique) + (1|census),
                           data=divers_woodsnow)

anova(mod.wsnow.shan1, mod.wsnow.shan2, mod.wsnow.shan3)
## refitting model(s) with ML (instead of REML)
## Data: divers_woodsnow
## Models:
## mod.wsnow.shan1: shannon ~ site + treatment + (1 | quad.unique) + (1 | census)
## mod.wsnow.shan3: shannon ~ treatment * season + (1 | quad.unique) + (1 | census)
## mod.wsnow.shan2: shannon ~ site + treatment + A.sum_curt + (1 | quad.unique) + 
## mod.wsnow.shan2:     (1 | census)
##                 Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## mod.wsnow.shan1  7 290.37 316.72 -138.18   276.37                        
## mod.wsnow.shan3  7 287.63 313.99 -136.82   273.63 2.734      0     <2e-16
## mod.wsnow.shan2  8 291.63 321.76 -137.82   275.63 0.000      1          1
##                    
## mod.wsnow.shan1    
## mod.wsnow.shan3 ***
## mod.wsnow.shan2    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod.wsnow.shan3, type=c('p', 'smooth'))

qqPlot(resid(mod.wsnow.shan3))

sjp.lmer(mod.wsnow.shan3, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wsnow.shan3, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.wsnow.shan3, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ treatment * season + (1 | quad.unique) + (1 | census)
##    Data: divers_woodsnow
## 
## REML criterion at convergence: 282.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.42405 -0.60104 -0.00628  0.59849  2.36740 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.18839  0.4340  
##  census      (Intercept) 0.11358  0.3370  
##  Residual                0.08021  0.2832  
## Number of obs: 319, groups:  quad.unique, 60; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          0.97980    0.21230   4.615
## treatmentS          -0.17670    0.12073  -1.464
## seasonsp            -0.22362    0.27856  -0.803
## treatmentS:seasonsp -0.25068    0.06445  -3.889
##----------------------- More diversity indices -------------------------##

### Invsimpson
range(divers_woodsnow$simpson)
## [1]   1 Inf
divers_woodsnow1 <- subset(divers_woodsnow, simpson != Inf)

mod.wsnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique)+(1|census),
                        data=divers_woodsnow1)

plot(mod.wsnow.simp, type=c('p', 'smooth'))

sjp.lmer(mod.wsnow.simp, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wsnow.simp, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.wsnow.simp))

summary(mod.wsnow.simp, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: divers_woodsnow1
## 
## REML criterion at convergence: 766.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50704 -0.54750 -0.00479  0.49057  2.62465 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.6946   0.8334  
##  census      (Intercept) 0.3738   0.6114  
##  Residual                0.4640   0.6812  
## Number of obs: 300, groups:  quad.unique, 59; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)           2.7789     0.3914   7.099
## site1                 0.4896     0.1633   2.998
## site2                 0.1122     0.1633   0.687
## treatmentS           -0.3765     0.2436  -1.545
## seasonsp             -0.4360     0.5103  -0.854
## treatmentS:seasonsp  -0.4229     0.1613  -2.621
### Pielou's evenness (J)
table(divers_woodsnow$species)
## 
##  0  1  2  3  4  5  6  7  8  9 
## 19 63 65 62 33 39 17 12  7  2
divers_woodsnow2 <- subset(divers_woodsnow, species > 1)
# 20% of the data don't meet the criteria
divers_woodsnow2$evenness <- divers_woodsnow2$shannon/log(divers_woodsnow2$species)
mod.wsnow.even <- lmer(evenness ~ site + treatment*season + (1|quad.unique) + (1|census),
                            data=divers_woodsnow2)

plot(mod.wsnow.even, type=c('p', 'smooth'))

sjp.lmer(mod.wsnow.even, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.wsnow.even, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

qqPlot(resid(mod.wsnow.even))

summary(mod.wsnow.even, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: evenness ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: divers_woodsnow2
## 
## REML criterion at convergence: -459.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6649 -0.6106  0.1042  0.6724  2.4046 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.004242 0.06513 
##  census      (Intercept) 0.003701 0.06083 
##  Residual                0.004832 0.06951 
## Number of obs: 237, groups:  quad.unique, 56; census, 6
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          0.885054   0.038181  23.180
## site1                0.032898   0.013948   2.359
## site2               -0.051215   0.013906  -3.683
## treatmentS           0.005687   0.021367   0.266
## seasonsp             0.029418   0.051147   0.575
## treatmentS:seasonsp -0.002994   0.019567  -0.153
###-------------------- Trees, shrubs -------------------------###
divers_snow_gf <- summarise(
  group_by(snowdat.ag_adult, site, quadrat, quad.unique, growth.form, treatment, census), 
  shannon=diversity(survs), simpson=diversity(survs, index='invsimpson'), species=specnumber(survs))

divers_snow_gf$A.sum <- adult.ngbr.inf$A.sum[match(divers_snow_gf$quad.unique, adult.ngbr.inf$quad.unique)]
divers_snow_gf$A.sum_curt <- divers_snow_gf$A.sum^(1/3)
summary(divers_snow_gf)
##  site       quadrat       quad.unique  growth.form treatment  census   
##  1:191   Min.   :102.0   1_1_305: 12   shrub:247   C:279     15sp: 86  
##  2:188   1st Qu.:203.0   1_1_306: 12   tree :266   S:234     15fa: 67  
##  3:134   Median :306.0   1_1_308: 12                         16sp: 56  
##          Mean   :299.9   1_1_401: 12                         16fa: 97  
##          3rd Qu.:404.0   1_1_501: 12                         17sp: 96  
##          Max.   :508.0   2_1_102: 12                         17fa:111  
##                          (Other):441                                   
##     shannon          simpson         species          A.sum     
##  Min.   :0.0000   Min.   :1.000   Min.   :0.000   Min.   :4487  
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:5335  
##  Median :0.5004   Median :1.800   Median :2.000   Median :5662  
##  Mean   :0.4698   Mean   :  Inf   Mean   :1.883   Mean   :5924  
##  3rd Qu.:0.7595   3rd Qu.:2.667   3rd Qu.:3.000   3rd Qu.:6533  
##  Max.   :1.7329   Max.   :  Inf   Max.   :6.000   Max.   :9223  
##                                                                 
##    A.sum_curt   
##  Min.   :16.49  
##  1st Qu.:17.47  
##  Median :17.82  
##  Mean   :18.05  
##  3rd Qu.:18.69  
##  Max.   :20.97  
## 
divers_snow_gf$season <- ifelse(divers_snow_gf$census=='15fa'| 
                                  divers_snow_gf$census=='16fa'|
                                  divers_snow_gf$census=='17fa', 'fa', 'sp')

summary(divers_snow_gf)
##  site       quadrat       quad.unique  growth.form treatment  census   
##  1:191   Min.   :102.0   1_1_305: 12   shrub:247   C:279     15sp: 86  
##  2:188   1st Qu.:203.0   1_1_306: 12   tree :266   S:234     15fa: 67  
##  3:134   Median :306.0   1_1_308: 12                         16sp: 56  
##          Mean   :299.9   1_1_401: 12                         16fa: 97  
##          3rd Qu.:404.0   1_1_501: 12                         17sp: 96  
##          Max.   :508.0   2_1_102: 12                         17fa:111  
##                          (Other):441                                   
##     shannon          simpson         species          A.sum     
##  Min.   :0.0000   Min.   :1.000   Min.   :0.000   Min.   :4487  
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:5335  
##  Median :0.5004   Median :1.800   Median :2.000   Median :5662  
##  Mean   :0.4698   Mean   :  Inf   Mean   :1.883   Mean   :5924  
##  3rd Qu.:0.7595   3rd Qu.:2.667   3rd Qu.:3.000   3rd Qu.:6533  
##  Max.   :1.7329   Max.   :  Inf   Max.   :6.000   Max.   :9223  
##                                                                 
##    A.sum_curt       season         
##  Min.   :16.49   Length:513        
##  1st Qu.:17.47   Class :character  
##  Median :17.82   Mode  :character  
##  Mean   :18.05                     
##  3rd Qu.:18.69                     
##  Max.   :20.97                     
## 
### Trees
# shannon diversity
mod.treesnow.shan <- lmer(shannon ~ site + treatment*season + (1|quad.unique) + (1|census),
                              data=subset(divers_snow_gf, growth.form=='tree'))

plot(mod.treesnow.shan, type=c('p', 'smooth'))

qqPlot(resid(mod.treesnow.shan))

sjp.lmer(mod.treesnow.shan, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.treesnow.shan, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.treesnow.shan, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: subset(divers_snow_gf, growth.form == "tree")
## 
## REML criterion at convergence: 206.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5386 -0.5750  0.0214  0.5469  2.6580 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.07709  0.2777  
##  census      (Intercept) 0.10313  0.3211  
##  Residual                0.07707  0.2776  
## Number of obs: 266, groups:  quad.unique, 60; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          0.57533    0.19489   2.952
## site1                0.12769    0.05645   2.262
## site2                0.02275    0.05633   0.404
## treatmentS          -0.01821    0.08564  -0.213
## seasonsp            -0.29592    0.26661  -1.110
## treatmentS:seasonsp -0.14373    0.07048  -2.039
## Invsimpson
divers_snow_gf1 <- subset(divers_snow_gf, simpson != Inf)
mod.treesnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique) + (1|census), 
                              data=subset(divers_snow_gf1, growth.form=='tree'))

plot(mod.treesnow.simp, type=c('p', 'smooth'))

qqPlot(resid(mod.treesnow.simp))

sjp.lmer(mod.treesnow.simp, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.treesnow.simp, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.treesnow.simp, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: subset(divers_snow_gf1, growth.form == "tree")
## 
## REML criterion at convergence: 462.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.92527 -0.55317 -0.08986  0.48331  3.11755 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.3505   0.5920  
##  census      (Intercept) 0.1581   0.3976  
##  Residual                0.2660   0.5158  
## Number of obs: 224, groups:  quad.unique, 59; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          1.90319    0.26199   7.264
## site1                0.38336    0.12017   3.190
## site2               -0.11325    0.11998  -0.944
## treatmentS           0.01019    0.17933   0.057
## seasonsp            -0.42786    0.33846  -1.264
## treatmentS:seasonsp -0.16242    0.15467  -1.050
### Shrubs
# shannon diversity
mod.shrubsnow.shan <- lmer(shannon ~ site + treatment*season + (1|quad.unique) + (1|census),
                               data=subset(divers_snow_gf, growth.form=='shrub'))

plot(mod.shrubsnow.shan, type=c('p', 'smooth'))

qqPlot(resid(mod.shrubsnow.shan))

sjp.lmer(mod.shrubsnow.shan, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.shrubsnow.shan, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.shrubsnow.shan, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: shannon ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: subset(divers_snow_gf, growth.form == "shrub")
## 
## REML criterion at convergence: 32.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.11544 -0.50524  0.01993  0.30473  3.13103 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.166842 0.40846 
##  census      (Intercept) 0.006976 0.08352 
##  Residual                0.030990 0.17604 
## Number of obs: 247, groups:  quad.unique, 52; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          0.54115    0.09828   5.506
## site1               -0.02548    0.08027  -0.317
## site2                0.06135    0.08153   0.752
## treatmentS          -0.36828    0.11839  -3.111
## seasonsp            -0.03984    0.07472  -0.533
## treatmentS:seasonsp -0.07002    0.04584  -1.527
## Invsimpson
mod.shrubsnow.simp <- lmer(simpson ~ site + treatment*season + (1|quad.unique) + (1|census), 
                               data=subset(divers_snow_gf1, growth.form=='shrub'))

plot(mod.shrubsnow.simp, type=c('p', 'smooth'))

qqPlot(resid(mod.shrubsnow.simp))

sjp.lmer(mod.shrubsnow.simp, type='fe')
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.

sjp.lmer(mod.shrubsnow.simp, type='re', sort.est="sort.all" )
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...
## Sorting each group of random intercept ('sort.all') is not possible when 'facet.grid = TRUE'.
## Plotting random effects...

summary(mod.shrubsnow.simp, correlation=FALSE)
## Linear mixed model fit by REML ['lmerMod']
## Formula: simpson ~ site + treatment * season + (1 | quad.unique) + (1 |  
##     census)
##    Data: subset(divers_snow_gf1, growth.form == "shrub")
## 
## REML criterion at convergence: 279.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9466 -0.3996 -0.0146  0.2641  3.8626 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  quad.unique (Intercept) 0.55828  0.7472  
##  census      (Intercept) 0.01539  0.1240  
##  Residual                0.08765  0.2961  
## Number of obs: 237, groups:  quad.unique, 51; census, 6
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)          1.88602    0.17157  10.993
## site1               -0.08559    0.14817  -0.578
## site2                0.12985    0.15021   0.864
## treatmentS          -0.64812    0.21748  -2.980
## seasonsp            -0.06676    0.11374  -0.587
## treatmentS:seasonsp -0.10162    0.07972  -1.275